Práctica Estadística Multivariante
Juan José Herrera Aranda, Pilar Navarro Ramírez, Adrián Rodríguez Montero
14 de enero de 2022
Análisis exploratorio univariante
Se ha realizado el análisis exploratorio univariante de la tercera base de datos. Primeramente, cargamos los datos del fichero DB_3.sav en un dataframe y con una codificación utf-8 y quitamos la primera columna que corresponde a los países, los cuales no nos dan información.
Observamos los datos con el objetivo de ver la forma en la que se nos presentan los datos, es decir, si son datos categóricos, numéricos continuos o discretos, positivos o negativos… y para ello se visualiza el dataframe. Observamos pues, que todos los datos son de tipo numérico continuo tanto positivos como negativos y cercanos entre sí. Destacar que las variables están estandarizadas. Además se puede apreciar que el atributo ZTLIBROP tiene un valor perdido en la fila 22.
## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legend## [1] "/home/pilarnavarro/pilar/sexto/Estadística Multivariante/prácticas/proyecto"datos3<-read.spss("DB_3.sav", to.data.frame=TRUE,reencode="utf-8")## re-encoding from utf-8datos3<-datos3[,-1]
datos3## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU
## 1 -0.857191432 0.96148387 -1.54985484 -0.21488857 -0.73382516 -0.606355478
## 2 -1.016717088 1.21338387 -0.89488793 -0.76526232 -0.88765979 0.057659570
## 3 -0.991855947 -0.31857942 0.43886285 1.13811356 1.26602514 -0.759937190
## 4 -1.077834060 -1.02030085 1.08192126 1.28029344 0.49685195 -1.058066395
## 5 -0.949384831 0.61190835 -0.32328046 0.44097348 -0.31077990 0.044108242
## 6 -1.072654656 -1.02801208 1.12955522 0.81706220 0.51608128 -1.107754596
## 7 -0.945241307 -0.46766309 0.47458832 1.06473039 -0.79151315 -0.606355478
## 8 -0.022271439 -0.29287533 0.27214400 -1.71924349 -0.80112781 1.905157220
## 9 2.861620943 -0.47280391 0.12924213 -0.04519000 -0.66652250 0.391925649
## 10 -0.667625230 1.48070632 -1.02588131 -0.60473665 -1.13764108 0.712640400
## 11 -0.326820419 -1.01001922 1.07001277 0.68405522 0.97758519 -0.497944858
## 12 0.571288308 0.04641855 -0.54954175 -0.93954734 -1.10879709 0.988184059
## 13 -0.052311985 -1.03829371 1.10573824 0.70240101 0.68914524 -0.958689993
## 14 0.094783101 -0.69899983 0.47458832 -0.20112923 1.46793310 -0.240469635
## 15 1.072654656 1.63493081 -1.45458692 -1.58165005 -0.86843046 1.801263710
## 16 -0.280205779 1.40359407 -1.54985484 -1.63668742 -1.14725575 1.232107954
## 17 -0.850976147 0.70958386 -0.58526722 -0.40293294 -1.07033843 0.292549247
## 18 0.035737891 -0.67072534 0.86756846 1.40871398 1.16987849 -1.071617722
## 19 0.842689099 -0.92776616 1.07001277 0.39510900 1.55446508 -0.705731880
## 20 2.152042537 -1.10255392 1.24864011 0.83999444 0.02573337 -0.881899137
## 21 1.561590433 -0.06153860 -0.06129369 0.80788931 -0.27232124 -0.705731880
## 22 -0.635512923 1.33162265 -1.01397283 -0.76067587 -1.10879709 0.708123290
## 23 -0.726670440 0.13895324 -0.07320218 0.39052255 -0.31077990 -0.186264325
## 24 -0.194227666 1.73774714 -2.14527929 -1.41195148 -1.13764108 1.728989963
## 25 0.001553821 1.90482367 -1.96665196 -1.36608700 -0.95496245 1.114663116
## 26 0.081316649 -0.74783758 0.65321565 0.01443382 0.92951186 -0.010097068
## 27 0.504991931 -0.98431514 0.78420903 0.83999444 1.27563980 -0.877382028
## 28 1.293297284 -0.98431514 0.98665335 1.50961584 0.01611870 -1.234233653
## 29 1.461109987 -0.96889269 0.91520242 1.21608317 0.93912653 -1.098720377
## 30 -0.128967170 -0.54220493 0.48649681 -0.44879742 0.48723728 -0.005579958
## 31 -0.508099574 1.06944101 -0.45427384 -0.64601468 -0.56076119 1.390206775
## 32 -0.974245972 -0.59361309 0.28405249 0.23916977 2.37171159 -0.439222439
## 33 -0.845796742 -0.97146310 0.97474486 0.72533325 0.68914524 -1.148408578
## 34 0.589934164 0.65817570 -0.79962002 -1.76969441 -1.00303577 1.832883474
## ZPSERVI ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## 1 0.48366865 -0.57510397 -0.56042975 -0.78289704 0.12813846
## 2 0.05831739 -0.96961561 -0.20336698 -2.13409403 -0.52613591
## 3 0.78468647 -0.47636769 -0.48303442 -0.41890863 -0.37922198
## 4 1.40635370 1.05949898 -0.61170369 0.62446672 1.18993245
## 5 0.11066832 -0.53207948 -0.71491390 -0.42726448 -0.72526157
## 6 1.62884513 1.57708364 -0.86586064 0.96600967 2.74978482
## 7 1.13151135 -0.69252676 -0.21044001 -0.71686695 -0.63345233
## 8 -1.88521069 -0.94691538 -0.38965874 1.70447158 -0.79477787
## 9 -0.50445505 0.77356425 1.16013447 -0.17435728 -0.51051973
## 10 -0.45864799 -0.88590734 0.24248176 -1.49621357 -0.82327121
## 11 0.26117722 1.42826571 0.10236959 -0.54331956 -0.16396932
## 12 -0.55680598 -0.93307433 -0.72408567 -0.43581567 -0.87196431
## 13 0.93519538 0.93292103 -0.02807325 0.37448801 0.55588195
## 14 -0.58298144 1.79320160 0.09042840 0.98956169 0.34057118
## 15 -1.57764901 -0.95544574 -0.73628387 -0.19176720 -0.92243683
## 16 -0.74657808 -0.96750778 -0.76279250 -0.33334671 -0.91277562
## 17 -0.67459556 -0.76440102 -0.07440193 -1.35279974 -0.70865028
## 18 1.35400278 0.99209789 4.42620176 -0.32337231 -0.19539657
## 19 0.43786159 -0.32373737 -0.25507169 -0.01157335 0.09841816
## 20 0.85666899 0.25235156 -0.74296346 0.94726825 0.31655099
## 21 1.13151135 -0.16469295 0.34209564 -1.43851975 -0.65883307
## 22 -0.76620967 NA -0.11195898 -1.25205613 -0.87754597
## 23 0.83703740 -0.94225917 -0.78996888 -0.98929097 -0.40072379
## 24 -1.36824531 -0.92401149 -0.69868856 -0.03897427 -0.94559840
## 25 -0.74657808 -0.96875479 -0.20915128 -1.20039046 -0.92516252
## 26 -0.63533236 0.07157343 -0.03681498 1.36222493 0.75427998
## 27 -0.16417404 0.44844158 0.04855887 1.65256659 1.56719986
## 28 1.21658160 1.75837557 -0.48919984 0.85527941 0.69960067
## 29 0.60145823 2.40236043 -0.10795229 0.46019499 1.05982477
## 30 -0.96252564 0.05946495 0.02611529 1.08311493 0.59174918
## 31 -1.17847321 -0.76872234 0.68707865 0.32878726 -0.74078430
## 32 0.01251033 -0.03092794 1.17613398 1.30446522 0.95562005
## 33 1.30819572 0.19806196 -0.42313576 0.88671706 2.65959005
## 34 -1.74778951 -0.92521141 1.92835266 0.72221179 -0.95066101En primer lugar mostramos un conjunto de estadísticos correspondientes a cada variable, la media, mediana, así como los valores máximos, mínimos, primer y tercer cuantil y existencia o no de valores perdidos. Aquí se puede apreciar de forma más rápida el único valor perdido presente en el conjunto de datos.
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB
## Min. :-1.0778 Min. :-1.1026 Min. :-2.1453 Min. :-1.7697
## 1st Qu.:-0.8497 1st Qu.:-0.9586 1st Qu.:-0.7460 1st Qu.:-0.7320
## Median :-0.1616 Median :-0.3931 Median : 0.2781 Median : 0.1268
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.5547 3rd Qu.: 0.8985 3rd Qu.: 0.9033 3rd Qu.: 0.8148
## Max. : 2.8616 Max. : 1.9048 Max. : 1.2486 Max. : 1.5096
##
## ZTMEDICO ZPAGRICU ZPSERVI ZTLIBROP
## Min. :-1.1473 Min. :-1.2342 Min. :-1.88521 Min. :-0.9696
## 1st Qu.:-0.8829 1st Qu.:-0.8480 1st Qu.:-0.72858 1st Qu.:-0.9240
## Median :-0.2916 Median :-0.2134 Median : 0.03541 Median :-0.3237
## Mean : 0.0000 Mean : 0.0000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.8694 3rd Qu.: 0.7115 3rd Qu.: 0.85176 3rd Qu.: 0.7736
## Max. : 2.3717 Max. : 1.9052 Max. : 1.62885 Max. : 2.4024
## NA's :1
## ZTEJERCI ZTPOBACT ZTENERGI
## Min. :-0.86586 Min. :-2.1341 Min. :-0.9507
## 1st Qu.:-0.59889 1st Qu.:-0.6735 1st Qu.:-0.7813
## Median :-0.20626 Median :-0.1067 Median :-0.3900
## Mean : 0.00000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.07996 3rd Qu.: 0.8789 3rd Qu.: 0.5828
## Max. : 4.42620 Max. : 1.7045 Max. : 2.7498
## Puesto que el dataset cuenta con valores perdidos, vemos en la siguiente tabla el número y la proporción de éstos con respecto al total de instancias del atributo respectivamente. Así pues, los valores perdidos de ZTLIBROP representan aproximadamente el 3% de la cantidad total, inferior al 5%, por lo que no se va a realizar ningún tipo de tratamiento a excepción de sustituirlo por la media correspondiente a dicho atributo.
cbind(apply(is.na(datos3),2,sum),apply(is.na(datos3),2,sum)/dim(datos3)[1])## [,1] [,2]
## ZPOBDENS 0 0.00000000
## ZTMINFAN 0 0.00000000
## ZESPVIDA 0 0.00000000
## ZPOBURB 0 0.00000000
## ZTMEDICO 0 0.00000000
## ZPAGRICU 0 0.00000000
## ZPSERVI 0 0.00000000
## ZTLIBROP 1 0.02941176
## ZTEJERCI 0 0.00000000
## ZTPOBACT 0 0.00000000
## ZTENERGI 0 0.00000000Comprobamos que efectivamente ha desaparecido el valor perdido y se ha sustituido correctamente por la media.
#Hay menos de 5% de los valores perdidos
not_available<-function(data,na.rm=F){
data[is.na(data)]<-mean(data,na.rm=T)
data
}
datos3$ZTLIBROP<-not_available(datos3$ZTLIBROP)Comprobamos si se han ido los valores perdidos
summary(datos3)## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB
## Min. :-1.0778 Min. :-1.1026 Min. :-2.1453 Min. :-1.7697
## 1st Qu.:-0.8497 1st Qu.:-0.9586 1st Qu.:-0.7460 1st Qu.:-0.7320
## Median :-0.1616 Median :-0.3931 Median : 0.2781 Median : 0.1268
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.5547 3rd Qu.: 0.8985 3rd Qu.: 0.9033 3rd Qu.: 0.8148
## Max. : 2.8616 Max. : 1.9048 Max. : 1.2486 Max. : 1.5096
## ZTMEDICO ZPAGRICU ZPSERVI ZTLIBROP
## Min. :-1.1473 Min. :-1.2342 Min. :-1.88521 Min. :-0.9696
## 1st Qu.:-0.8829 1st Qu.:-0.8480 1st Qu.:-0.72858 1st Qu.:-0.9145
## Median :-0.2916 Median :-0.2134 Median : 0.03541 Median :-0.2442
## Mean : 0.0000 Mean : 0.0000 Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.8694 3rd Qu.: 0.7115 3rd Qu.: 0.85176 3rd Qu.: 0.6923
## Max. : 2.3717 Max. : 1.9052 Max. : 1.62885 Max. : 2.4024
## ZTEJERCI ZTPOBACT ZTENERGI
## Min. :-0.86586 Min. :-2.1341 Min. :-0.9507
## 1st Qu.:-0.59889 1st Qu.:-0.6735 1st Qu.:-0.7813
## Median :-0.20626 Median :-0.1067 Median :-0.3900
## Mean : 0.00000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.07996 3rd Qu.: 0.8789 3rd Qu.: 0.5828
## Max. : 4.42620 Max. : 1.7045 Max. : 2.7498Exploración descriptiva
En este apartado iremos variable por variable obteniendo los resultados de aplicar diferentes medidas descriptivas, clásicas y resistentes, de centralidad, forma y dispersión.
Variable ZPOBDENS
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 1.06373e-16
##
##
## $centralidad$resistente
## [,1]
## PMC -0.1474835
## Trimedia -0.1545405
## Centrimedia -0.2293829
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 9.400883e+15
##
## $dispersion$clasica$rango
## [1] -1.077834 2.861621
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.4043955
## MEDA 0.6924864
## CVc -4.7611940
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 1.067954
##
## $forma$clasica$kurtosis
## [1] 0.5421592
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.08733974
## Asimetría de Kelly -0.37369403
## Asimetría de Kelly adimensional 2.31250000
Las medidas resistentes de centralidad están ligeramente desplazadas hacia la izquierda siendo PMC muy próxima a cero. El valor de MEDA es ligeramente inferior al de la desviación típica. En cuanto a las medidas de simetría tenemos que el coeficiente skewness es positivo, por lo que la distribución tendrá una cola asimétrica extendida hacia los valores positivos. Tenemos una curtosis positiva, por lo que la distribución es leptocúrtica habiendo una mayor concentración de datos entorno a la media.
En el histograma vemos una forma bastante anormal para ser una distribución estandarizada.
Variable ZTMINFAN
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 8.982296e-17
##
##
## $centralidad$resistente
## [,1]
## PMC -0.0300511
## Trimedia -0.2115862
## Centrimedia -0.2268480
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 1.113301e+16
##
## $dispersion$clasica$rango
## [1] -1.102554 1.904824
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.8571199
## MEDA 0.6040459
## CVc -30.8993711
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.5477565
##
## $forma$clasica$kurtosis
## [1] -1.203999
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.9235577
## Asimetría de Kelly -0.6132994
## Asimetría de Kelly adimensional 1.5600769
Las medidas resistentes de centralidad están ligeramente desplazadas hacia la izquierda siendo PMC muy próxima a cero. El valor de MEDA es inferior al de la desviación típica. En cuanto a las medidas de simetría tenemos que el coeficiente skewness es positivo, por lo que la distribución tendrá una cola asimétrica extendida hacia los valores positivos. Tenemos una curtosis negativa, por lo que la distribución es platicúrtica habiendo una menor concentración de datos entorno a la media.
En el histograma se puede observar que lo anterior, es debido a que se tienen dos subgrupos de datos.
Variable ZESPVIDA
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] -1.72884e-16
##
##
## $centralidad$resistente
## [,1]
## PMC 0.07863105
## Trimedia 0.17836465
## Centrimedia 0.17613181
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] -5.784225e+15
##
## $dispersion$clasica$rango
## [1] -2.145279 1.248640
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.6493257
## MEDA 0.7621433
## CVc 10.4877506
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.5802673
##
## $forma$clasica$kurtosis
## [1] -0.843727
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.7172544
## Asimetría de Kelly 0.4995611
## Asimetría de Kelly adimensional 1.7963476
Las medidas resistentes de centralidad están suavemente desplazadas a la derecha aproximandose el PMC a cero. Lo anterior sumado al rango indica un amontonamiento de los datos en la parte que está a la derecha del cero. El MEDA está ligeramente por debajo a la desviación típica. En cuanto a los estimadores de asimetría tenemos un coeficiente Skewness negativo, por lo que la distribución tendrá una cola asimétrica extendida hacia los valores negativos. Además, tenemos un coeficiente de curtosis también negativo por lo que la distribición es platicúrtica, presentando así una menor distribución de datos en torno a la media. Lo anterior se puede ver de manera visual en el histograma.
Variable ZPOBURB
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 4.511206e-16
##
##
## $centralidad$resistente
## [,1]
## PMC 0.0413792
## Trimedia 0.0840905
## Centrimedia 0.1147624
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 2.216702e+15
##
## $dispersion$clasica$rango
## [1] -1.769694 1.509616
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.5467796
## MEDA 0.7223655
## CVc 18.6903015
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.318985
##
## $forma$clasica$kurtosis
## [1] -1.097519
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.6736702
## Asimetría de Kelly 0.2958259
## Asimetría de Kelly adimensional 2.3329787
Las medidas de centralidad están un poco desplazadas ligeramente hacia la derecha que junto con el rango tenemos que los datos se agrupan más en una región centrada ligeramente desviada hacia la derecha del cero. EL MEDA se queda por debajo de la desviación típica. En cuanto a los coeficientes de asimetría tenemos un skewness negativo por lo que la distribución tendrá una cola asimétrica extendida hacia los valores negativos y la curtosis es menor que cero, así pues estamos ante una distribución platicúrtica. Lo dicho anteriormente se observa de manera intuitiva en la gráfica.
El histograma parece acorde con una normal.
Variable ZTMEDICO
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 2.416495e-17
##
##
## $centralidad$resistente
## [,1]
## PMC -0.006716126
## Trimedia -0.149133349
## Centrimedia -0.149734266
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 4.138224e+16
##
## $dispersion$clasica$rango
## [1] -1.147256 2.371712
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.7522727
## MEDA 0.7980172
## CVc -130.4526316
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.4909934
##
## $forma$clasica$kurtosis
## [1] -0.9345197
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.9769641
## Asimetría de Kelly -0.3735297
## Asimetría de Kelly adimensional 1.2811833
Las medidas resistentes de centralidad están ligeramente desplazadas hacia la izquierda que sumado al rango indican que los datos se distribuyen en una amplia zona situados a la izquierda del cero. El MEDA se queda por debajo de la desviación típica. En cuento a los coeficientes de simetría tenemos que el índice skewness es positivo por lo que la distribución tendrá una cola asimétrica extendida hacia los valores positivos. Además la curtosis es negativa, por lo que la disribución es platicúrtica.
Variable ZPAGRICU
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 6.355338e-17
##
##
## $centralidad$resistente
## [,1]
## PMC -0.06825485
## Trimedia -0.14081091
## Centrimedia -0.20433276
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 1.57348e+16
##
## $dispersion$clasica$rango
## [1] -1.234234 1.905157
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.5595319
## MEDA 0.7069276
## CVc -11.4243309
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.5918253
##
## $forma$clasica$kurtosis
## [1] -0.9443355
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.6801059
## Asimetría de Kelly -0.4817497
## Asimetría de Kelly adimensional 2.2578456
Las medidas resistentes de centralidad están ligeramente desplazadas hacia la izquierda siendo PMC muy próxima a cero. El valor de MEDA es ligeramente inferior al de la desviación típica. En cuanto a las medidas de simetría tenemos que el coeficiente skewness es positivo, por lo que la distribución tendrá una cola asimétrica extendida hacia los valores positivos. Tenemos una curtosis negativa, por lo que la distribución es platicúrtica habiendo una menor concentración de datos entorno a la media. Todo lo anterior se puede apreciar en el histograma.
Variable ZPSERVI
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] -3.529078e-17
##
##
## $centralidad$resistente
## [,1]
## PMC 0.0615893232
## Trimedia 0.0485015920
## Centrimedia 0.0006495749
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] -2.833601e+16
##
## $dispersion$clasica$rango
## [1] -1.885211 1.628845
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.5803435
## MEDA 0.7918077
## CVc 12.8296875
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.1239223
##
## $forma$clasica$kurtosis
## [1] -1.088712
##
##
## $forma$resistente
## 25%
## Asimetría de Yule 0.73913043
## Asimetría de Kelly 0.05071496
## Asimetría de Kelly adimensional 1.43206522
Las medidas de centralidad están un poco desplazadas hacia la derecha que junto con el rango tenemos que los datos se agrupan más en una región centrada ligeramente desviada hacia la derecha del cero. EL MEDA se queda por debajo de la desviación típica. En cuanto a los coeficientes de asimetría tenemos un skewness negativo por lo que la distribución tendrá una cola asimétrica extendida hacia los valores negativos y la curtosis es menor que cero, así pues estamos ante una distribución platicúrtica. Lo dicho anteriormente se observa de manera intuitiva en la gráfica.
Variable ZTLIBROP
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 6.729069e-17
##
##
## $centralidad$resistente
## [,1]
## PMC -0.1111009
## Trimedia -0.1776580
## Centrimedia -0.2615358
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 0.9847319
##
## $dispersion$clasica$Coef_varización
## [1] 1.4634e+16
##
## $dispersion$clasica$rango
## [1] -0.9696156 2.4023604
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.6067690
## MEDA 0.6849277
## CVc -7.2311230
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 0.8170779
##
## $forma$clasica$kurtosis
## [1] -0.4758477
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.5450695
## Asimetría de Kelly -0.5339910
## Asimetría de Kelly adimensional 2.1865595
Las medidas resistentes de centralidad están ligeramente desplazadas hacia la izquierda siendo PMC próxima a cero. El valor de MEDA es inferior al de la desviación típica. En cuanto a las medidas de simetría tenemos que el coeficiente skewness es positivo, por lo que la distribución tendrá una cola asimétrica extendida hacia los valores positivos. Tenemos una curtosis negativa, por lo que la distribución es platicúrtica habiendo una menor concentración de datos entorno a la media.
En el histograma vemos una forma bastante anormal para ser una distribución estandarizada.
Variable ZTEJERCI
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 1.859184e-16
##
##
## $centralidad$resistente
## [,1]
## PMC -0.2594621
## Trimedia -0.2328606
## Centrimedia -0.2192510
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 5.378705e+15
##
## $dispersion$clasica$rango
## [1] -0.8658606 4.4262018
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 0.6788462
## MEDA 0.3313997
## CVc -1.3081800
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 2.86539
##
## $forma$clasica$kurtosis
## [1] 9.634918
##
##
## $forma$resistente
## 25%
## Asimetría de Yule 0.2579423
## Asimetría de Kelly -0.3448882
## Asimetría de Kelly adimensional 1.6721112
Las medidas resistentes de centralidad están desplazadas hacia la izquierda siendo PMC muy próxima a cero. El valor de MEDA es ligeramente inferior al de la desviación típica. En cuanto a las medidas de simetría tenemos que el coeficiente skewness es positivo, por lo que la distribución tendrá una cola asimétrica extendida hacia los valores positivos. Tenemos una curtosis positiva, por lo que la distribución es leptocúrtica habiendo una mayor concentración de datos entorno a la media.
En el histograma vemos una forma bastante anormal para ser una distribución estandarizada. Además podemos observar que muy probablemente existan outliers.
Variable ZTPOBACT
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] -2.952027e-16
##
##
## $centralidad$resistente
## [,1]
## PMC 0.10268877
## Trimedia -0.00198850
## Centrimedia 0.02917055
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] -3.387502e+15
##
## $dispersion$clasica$rango
## [1] -2.134094 1.704472
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.5523378
## MEDA 0.8557514
## CVc 7.5584589
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] -0.1284517
##
## $forma$clasica$kurtosis
## [1] -0.8884493
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -1.96271531
## Asimetría de Kelly -0.06440751
## Asimetría de Kelly adimensional 0.60382547
Las medidas de centralidad están muy próximas a cero, siendo PMC positiva y muy próxima a cero. EL valor de MEDA se queda por debajo de la desviación típica. En cuanto a los coeficientes de asimetría tenemos un skewness negativo próximo a cero por lo que la distribución tendrá una cola asimétrica extendida hacia los valores negativos y la curtosis es menor que cero, así pues estamos ante una distribución platicúrtica. Lo dicho anteriormente se observa de manera intuitiva en la gráfica.
El histograma parece acorde con una normal.
Variable ZTENERGI
## $centralidad
## $centralidad$clasica
## $centralidad$clasica$media
## [1] 3.586955e-17
##
##
## $centralidad$resistente
## [,1]
## PMC -0.09924855
## Trimedia -0.24461072
## Centrimedia -0.26271176
##
##
## $dispersion
## $dispersion$clasica
## $dispersion$clasica$desviación_típica
## [1] 1
##
## $dispersion$clasica$Coef_varización
## [1] 2.78788e+16
##
## $dispersion$clasica$rango
## [1] -0.950661 2.749785
##
##
## $dispersion$resistente
## 75%
## Rango Inter-Cuartílico 1.364062
## MEDA 0.520457
## CVc -6.871948
##
##
## $forma
## $forma$clasica
## $forma$clasica$skewness
## [1] 1.227034
##
## $forma$clasica$kurtosis
## [1] 0.9124716
##
##
## $forma$resistente
## 25%
## Asimetría de Yule -0.7454988
## Asimetría de Kelly -0.5056537
## Asimetría de Kelly adimensional 1.2966382
Las medidas resistentes de centralidad están ligeramente desplazadas hacia la izquierda siendo PMC muy próxima a cero. El valor de MEDA es inferior al de la desviación típica. En cuanto a las medidas de simetría tenemos que el coeficiente skewness es positivo, por lo que la distribución tendrá una cola asimétrica extendida hacia los valores positivos. Tenemos una curtosis positiva, por lo que la distribución es leptocúrtica habiendo una mayor concentración de datos entorno a la media.
El histograma es de una distribución que se aleja bastante de una normal. De nuevo tenemos la sospecha de que existan outliers en esta variable.
Por último, en este apartado presentamos de forma conjunta los histogramas de todas las variables:
Exploración Gráfica
Una forma de comprobar si hay outliers es representar un diagrama de cajas el cual para el conjunto de datos escogido es el siguiente.
boxplot(datos3,main="Datos originales" ,
xlab= "",
ylab= "z-values",
col=c(1:11),
las=2)
Podemos comprobar que sí hay outliers. Puesto que no disponemos de demasiados datos en nuestra muestra (únicamente 34), no es conveniente eliminar los outliers, de modo que los sustituiremos por la media. Ahora definimos la función outliers para sustituir estos valores por la media. Debido a la gran variabilidad de los atributos, ejecutaremos la función varias veces para poder eliminar los outliers.
outlier<-function(data,na.rm=T){
H<-1.5*IQR(data)
if(any(data<=(quantile(data,0.25,na.rm = T)-H))){
data[data<=quantile(data,0.25,na.rm = T)-H]<-NA
data[is.na(data)]<-mean(data,na.rm=T)
data<-outlier(data)}
if(any(data>=(quantile(data,0.75, na.rm = T)+H))){
data[data>=quantile(data,0.75, na.rm = T)+H]<-NA
data[is.na(data)]<-mean(data,na.rm=T)
data<-outlier(data)
}
return(data)
}El resultado final es el siguiente, donde podemos ver que ya no hay outliers:
Presentamos ambos gráficos juntos para comparar los datos originales y los arreglados.
par(mfrow=c(1,2))
# Boxplot de los datos originales
boxplot(datos3, main="Datos originales",
xlab="Variables",
ylab="z-values",
col=c(1:11))
# Boxplot de los datos corregidos.
boxplot(datos3_nout,main="Datos sin outliers",
xlab="Variables",
ylab="z-values",
col=c(1:11))
Normalidad
Para poder aplicar ciertas técnicas estadísticas, necesitamos saber si estamos tratando con variables normales, para ello usaremos el método gráfico de comparación de cuantiles qq-plot.
Datos sin outliers
par(mar=c(1,1,1,1))
par(mfrow=c(3,5))
invisible(apply(datos3_nout, 2, function(x){
qqnorm(x,main=NULL)
abline(a=0,b=1,col="red")
}))
Datos originales
par(mar=c(1,1,1,1))
par(mfrow=c(3,5))
invisible(apply(datos3, 2, function(x){
qqnorm(x,main=NULL)
abline(a=0,b=1,col="red")
}))
Atendiendo a los gráficos presentados vemos como las variables que mas se acercan a la normalidad son la 4, la 6, la 7 y la 10. No tan cerca están las variables 1, 2, 3, 5, 8, 11, y con desviaciones totales encontramos la variable 9 (estas desviaciones se pueden deber a los valores outliers no tratados). Para estar más seguros vamos a realizar el test de Kolmogorov-Smirnov para la variable 9 y para algunas de las variables que presentan una pequeña desviación.
Test de Kolmogorov-Smirnov:
Vamos a aplicar el test de Kolmogorov-Smirnov para probar la bondad del ajuste, es decir, para ver si la distribución muestral se ajusta a una distribución normal. La hipótesis nula establece que la distribución empírica es similar a una distribución normal frente a la alternativa que establece que la distribución de frecuencias observada no es consistente con la distribución normal.
Para el caso de la variable ZTEJERCI vemos que el p-valor obtenido es 0.02069, por tanto al 95% de confianza rechazamos la hipótesis nula, es decir, no seguiría una distribución normal. Si aumentamos el nivel de confianza hasta el 99% no podemos rechazar la hipótesis nula y, por tanto se podría admitir la normalidad.
ks.test(datos3$ZTEJERCI,"pnorm",0,1)##
## One-sample Kolmogorov-Smirnov test
##
## data: datos3$ZTEJERCI
## D = 0.25335, p-value = 0.02069
## alternative hypothesis: two-sidedHemos realizado el test sobre algunas de las variables que presentan una ligera desviación atendiendo al método gráfico qq-plot y hemos podido comprobar que en el resto de variables no podemos rechazar la hipótesis nula ya que, el p-valor es superior al valor de 0.2 en todas las variables.
ks.test(datos3$ZPOBDENS,"pnorm",0,1)##
## One-sample Kolmogorov-Smirnov test
##
## data: datos3$ZPOBDENS
## D = 0.16813, p-value = 0.261
## alternative hypothesis: two-sidedks.test(datos3$ZTLIBROP,"pnorm",0,1)##
## One-sample Kolmogorov-Smirnov test
##
## data: datos3$ZTLIBROP
## D = 0.16612, p-value = 0.2736
## alternative hypothesis: two-sidedks.test(datos3$ZTENERGI,"pnorm",0,1)##
## One-sample Kolmogorov-Smirnov test
##
## data: datos3$ZTENERGI
## D = 0.17715, p-value = 0.2098
## alternative hypothesis: two-sidedCorrelación entre las variables
##
## Attaching package: 'psych'## The following object is masked _by_ '.GlobalEnv':
##
## outlier## Loading required package: ggplot2##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alpha## corrplot 0.90 loaded##
## Attaching package: 'polycor'## The following object is masked from 'package:psych':
##
## polyserialEn primer lugar vemos la matriz de correlación, para tener una idea de la correlación en la muestra.
cor(datos_pca)## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU
## ZPOBDENS 1.00000000 -0.2109635 0.1735014 0.05132529 0.0525974 -0.018865672
## ZTMINFAN -0.21096354 1.0000000 -0.9668834 -0.75749544 -0.7509189 0.752161039
## ZESPVIDA 0.17350137 -0.9668834 1.0000000 0.78714604 0.7361282 -0.753537421
## ZPOBURB 0.05132529 -0.7574954 0.7871460 1.00000000 0.6353621 -0.938022359
## ZTMEDICO 0.05259740 -0.7509189 0.7361282 0.63536205 1.0000000 -0.674513147
## ZPAGRICU -0.01886567 0.7521610 -0.7535374 -0.93802236 -0.6745131 1.000000000
## ZPSERVI -0.05804281 -0.5901653 0.6121270 0.89001519 0.4446255 -0.907218449
## ZTLIBROP 0.23071577 -0.7157201 0.7004692 0.66106622 0.6093143 -0.666850181
## ZTEJERCI 0.05154113 -0.0254503 0.1017843 0.03566091 0.1513290 -0.007128371
## ZTPOBACT 0.23749261 -0.6032470 0.5411995 0.15488847 0.5336484 -0.147150671
## ZTENERGI 0.17123880 -0.7155678 0.6651297 0.61549363 0.7507638 -0.688730826
## ZPSERVI ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## ZPOBDENS -0.05804281 0.2307158 0.051541133 0.23749261 0.1712388
## ZTMINFAN -0.59016525 -0.7157201 -0.025450301 -0.60324698 -0.7155678
## ZESPVIDA 0.61212697 0.7004692 0.101784303 0.54119951 0.6651297
## ZPOBURB 0.89001519 0.6610662 0.035660910 0.15488847 0.6154936
## ZTMEDICO 0.44462549 0.6093143 0.151329048 0.53364837 0.7507638
## ZPAGRICU -0.90721845 -0.6668502 -0.007128371 -0.14715067 -0.6887308
## ZPSERVI 1.00000000 0.5044833 -0.173017490 -0.05616141 0.4159030
## ZTLIBROP 0.50448334 1.0000000 0.113918652 0.41591469 0.6602394
## ZTEJERCI -0.17301749 0.1139187 1.000000000 -0.07216152 0.1076248
## ZTPOBACT -0.05616141 0.4159147 -0.072161520 1.00000000 0.6004197
## ZTENERGI 0.41590296 0.6602394 0.107624821 0.60041968 1.0000000Representamos visualmente las correlaciones:
poly_cor<-hetcor(datos_pca)$correlations
ggcorrplot(poly_cor, type="lower",hc.order=T)
corrplot(cor(datos_pca), order = "hclust", tl.col='black', tl.cex=1)
correlaciones <- correlate(datos_pca) ##
## Correlation method: 'pearson'
## Missing treated using: 'pairwise.complete.obs'rplot(correlaciones, legend = TRUE, colours = c("firebrick1", "black","darkcyan"), print_cor = TRUE) ## Don't know how to automatically pick scale for object of type noquote. Defaulting to continuous.## Don't know how to automatically pick scale for object of type noquote. Defaulting to continuous.
Podemos observar que hay cierta correlación entre las variables, ya que la matriz de correlaciones dista de la matriz identidad.
Por ejemplo vemos que la variable ZTMINFAN esta correlada con la variable ZESPVIDA en un 96.7 % o la variable ZPAGRICU con ZPOBURB en un 93.8%.
Para comprobar que a nivel poblacional efectivamente están correladas las variables usaremos el test de Bartlett. La hipótesis nula es que el determinante de la matriz de correlaciones es uno.
Para aplicar el test de Bartlett, primero es necesario normalizar los datos:
# Se normalizan los datos
datos_normalizados<-scale(datos_pca)
# Se hace el test de esfericidad
cortest.bartlett(cor(datos_normalizados),n=ncol(datos_normalizados))## $chisq
## [1] 74.79283
##
## $p.value
## [1] 0.03913137
##
## $df
## [1] 55Como el p-valor es 0.03913137 (inferior 0,15) rechazamos la hipótesis nula y concluimos que las variables están correladas, luego tiene sentido aplicar análisis de componentes principales.
Análisis de componentes principales (ACP)
Llevamos a cabo un análisis de componentes principales y se muestran la desviación típica, la varianza explicada y varianza acumulada de dichas componentes.
# Pasamos los par?metros "scale" y "center" a TRUE para considerar
# los datos originales normalizados
PCA<-prcomp(datos_pca, scale=T, center = T)
#Obtenemos la desviación típica, la varianza explicada y varianza acumulada de las componentes principales
summary(PCA)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.4818 1.2806 1.03112 0.95154 0.6448 0.60272 0.48187
## Proportion of Variance 0.5599 0.1491 0.09665 0.08231 0.0378 0.03302 0.02111
## Cumulative Proportion 0.5599 0.7090 0.80568 0.88799 0.9258 0.95881 0.97992
## PC8 PC9 PC10 PC11
## Standard deviation 0.33542 0.24916 0.16420 0.13908
## Proportion of Variance 0.01023 0.00564 0.00245 0.00176
## Cumulative Proportion 0.99015 0.99579 0.99824 1.00000Vemos la proporción de varianza explicada:
varianza_explicada <- PCA$sdev^2 / sum(PCA$sdev^2)
varianza_explicada## [1] 0.559940548 0.149080884 0.096654411 0.082312012 0.037797654 0.033024464
## [7] 0.021108960 0.010227819 0.005643820 0.002450981 0.001758449#Gráficamente
ggplot(data = data.frame(varianza_explicada, pc = 1:11),
aes(x = pc, y = varianza_explicada, fill=varianza_explicada )) +
geom_col(width = 0.3) +
scale_y_continuous(limits = c(0,0.6)) + theme_bw() +
labs(x = "Componente principal", y= " Proporción de varianza explicada")
Vemos la varianza acumulada:
varianza_acum<-cumsum(varianza_explicada)
varianza_acum## [1] 0.5599405 0.7090214 0.8056758 0.8879879 0.9257855 0.9588100 0.9799189
## [8] 0.9901468 0.9957906 0.9982416 1.0000000#Gráficamente
ggplot( data = data.frame(varianza_acum, pc = 1:11),
aes(x = pc, y = varianza_acum ,fill=varianza_acum )) +
geom_col(width = 0.5) +
scale_y_continuous(limits = c(0,1)) +
theme_bw() +
labs(x = "Componente principal",
y = "Proporción varianza acumulada")
Selección del número óptimo de componentes principales
Usamos la regla de Abdi et al. (2010), esto es, se promedian las varianzas de las componentes principales y se seleccionan aquellas cuya proporción de varianza explicada supera la media.
## [1] 6.15934602 1.63988972 1.06319852 0.90543213 0.41577420 0.36326910
## [7] 0.23219856 0.11250601 0.06208202 0.02696079 0.01934294## [1] 1Como la media es 1, nos quedamos con las tres primeras componentes principales, pues la proporción de varianza explicada por cada una de ellas es mayor a 1. La varianza acumulada es 0.8056758, luego estas componentes pricipales explican el 80,57% de la varianza original.
Además, usamos ahora el método del codo con Scree plot:
scree(poly_cor) 
En efecto, el método del codo confirma lo que ya sabíamos, que el número óptimo de componentes principales es 3.
Pesos de cada variable en el sistema determinado por las tres primeras componentes principales
Veamos ahora los pesos de cada variable en la correspondiente componente principal.
## PC1 PC2 PC3
## ZPOBDENS 0.06796131 0.39036079 0.07578571
## ZTMINFAN -0.37529798 -0.11573027 -0.08212636
## ZESPVIDA 0.37165513 0.07270402 -0.00407411
## ZPOBURB 0.35808785 -0.29124050 -0.06094609
## ZTMEDICO 0.33374417 0.13909259 -0.09587119
## ZPAGRICU -0.36224935 0.29940891 0.04348437
## ZPSERVI 0.29607632 -0.48878963 0.09183567
## ZTLIBROP 0.32594252 0.08892433 -0.07078859
## ZTEJERCI 0.02397325 0.17948476 -0.93319975
## ZTPOBACT 0.20289865 0.56167354 0.29659750
## ZTENERGI 0.33155119 0.20152127 -0.02054058Así, si consideramos la base formada por las variables {ZPOBDENS, ZTMINFAN, ZESPVIDA, ZPOBURB, ZTMEDICO, ZPAGRICU, ZPSERVI, ZTLIBROP, ZTEJERCI, ZTPOBACT, ZTENERGI} las coordenadas de cada una de las componentes principales en dicha base serían:
## [1] "PC1"## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU
## 0.06796131 -0.37529798 0.37165513 0.35808785 0.33374417 -0.36224935
## ZPSERVI ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## 0.29607632 0.32594252 0.02397325 0.20289865 0.33155119## [1] "PC2"## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU
## 0.39036079 -0.11573027 0.07270402 -0.29124050 0.13909259 0.29940891
## ZPSERVI ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## -0.48878963 0.08892433 0.17948476 0.56167354 0.20152127## [1] "PC3"## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU
## 0.07578571 -0.08212636 -0.00407411 -0.06094609 -0.09587119 0.04348437
## ZPSERVI ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## 0.09183567 -0.07078859 -0.93319975 0.29659750 -0.02054058Las coordenadas de los 6 primeros datos originales (tipificados) en el sistema de referencia determinado por las 3 componentes principales son las siguientes:
## PC1 PC2 PC3
## [1,] -1.1879510 -1.5717487 0.4504449
## [2,] -2.3462462 -1.8810662 -0.8108192
## [3,] 1.1979236 -1.5383367 0.2007784
## [4,] 3.2068509 -0.8941585 0.8394015
## [5,] -0.8909634 -1.3473518 0.8989330
## [6,] 2.7566836 -1.1358451 1.6338045## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBaPodríamos mostrar las coordenadas de las 34 observaciones, pero con las 6 primeras es suficiente para hacernos una idea del susodicho sistema de referencia.
Ahora veamos gráficamente el peso que tienen las variables en la definición de las componentes principales.
Comparativa entre la primera y segunda componente principal:
fviz_pca_var(PCA,
repel=TRUE,col.var="cos2",
legend.title="Distancia")+theme_bw()
Vemos que algunas variables como ZPSERVI o ZTPOBACT están casi correladas por igual con ambas componentes principales. Sin embargo, hay muchas que están más correladas con la primera componente principal como ZTLIBROP, ZTMEDICO y ZESPVIDA, siendo esta última la que presenta una mayor correlación positiva con dicha componente principal, mientras que ZTMINFAN presenta una mayor correlación negativa con la primera componente. Por otro lado tenemos un par de variables, ZTEJERC y ZPOBDENS, que están más correladas con la segunda componente principal.
Comparativa entre la primera y tercera componente principal:
fviz_pca_var(PCA,axes=c(1,3),
repel=TRUE,col.var="cos2",
legend.title="Distancia")+theme_bw()
En este caso la mayoría de las variables están más correladas con la primera componente principal, tanto positiva como negativamente. únicamente vemos una variable, ZTEJERCI que está casi perfectamente correlada negativamente con la tercera componente principal. Apreciamos además una variable, ZTPOBACT, que presenta prácticamente la misma correlación positiva con ambas componentes principales. Además, notamos que hay una mayor correlación de las variables con la primera componente principal respecto a la tercera de la que había con la primera componente con respecto a la segunda en el primer gráfico, ya que las flechas están más pegadas al eje horizontal.
Comparativa entre la segunda y tercera componente principal
fviz_pca_var(PCA,axes=c(2,3),
repel=TRUE,col.var="cos2",
legend.title="Distancia")+theme_bw()
Como ocurría en el caso anterior, solo la variable ZTEJERCI está fuertemente correlada con la tercera componente principal, aunque aquí esta variable tiene algo más de peso en la segunda componente principal de la que tenía en la primera. Vemos que el resto de variables están más correladas con la segunda componente principal, aunque en general no en gran medida, pues las flechas son pequeñas y de color oscuro. La variable ZTPOBACT, al igual que en el gráfico anterior, influye casi por igual en ambas componentes principales, siendo en este gráfico la correlación positiva mayor, pues la flecha es de mayor longitud.
Podemos concluir entonces que la componente principal con la que están menos correladas las variables es la tercera, siendo únicamente ZTEJERCI la variable que más peso tiene en la tercera componente, aunque el resto también influyen un poco en dicha componente. Por ejemplo ZTPOBACT tiene un peso considerable en la tercera componente. El resto de variables tienen, en general, más peso en la primera componente, aunque su influencia en la segunda componente también es considerable, estando algunas variables más correladas con la segunda que con la primera.
Representamos ahora conjuntamente las variables y las observaciones relacionando visualmente las posibles relaciones entre las observaciones, las contribuciones de estas a las varianzas de las componentes y el peso de las variables en cada componentes principal. Con la orden “contrib” podemos identificar con colores aquellas observaciones que mayor varianza explican de las componentes principales, siendo la influencia de la observación menor cuanto más cerca se encuentre esta del origen.
Variables y observaciones en la 1º y 2ª componente principal
fviz_pca(PCA,axes=c(1,2),
alpha.ind ="contrib", col.var = "cos2",col.ind="seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
repel=TRUE,
legend.title="Distancia")+theme_bw()
Notamos que la observación que menos influye en la varianza explicada es la 9 pues está prácticamente sobre el origen de coordenadas. Las observaciones que más afectan a la primera componente principal son la 24 y la 29, ya que están más en los extremos del eje horizontal. Por otra parte, en la segunda componente principal, tienen más peso la observación 7 y 8.
Variables y observaciones en la 1º y 3ª componente principal
fviz_pca(PCA,axes=c(1,3),
alpha.ind ="contrib", col.var = "cos2",col.ind="seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
repel=TRUE,
legend.title="Distancia")+theme_bw()
En este caso observamos que los datos 9 y 7 son los que menos afectan a las primera y tercera componente principal. Al igual que observamos en el gráfico anterior, los datos 29 y 24 son los que mayor influencia tienen en la primera componente, mientras que las observaciones 20 y 31 afectan en mayor medida a la tercera componente.
Variables y observaciones en la 2ª y 3ª componente principal
fviz_pca(PCA,axes=c(2,3),
alpha.ind ="contrib", col.var = "cos2",col.ind="seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"),
repel=TRUE,
legend.title="Distancia")+theme_bw()
La observación 9 vuelve a ser la que menos afecta a las componentes principales, seguida de la 19 y la 13. La segunda componente se ve más influenciada, como ya sabíamos, por las observaciones 7 y 8. Con respecto a la tercera componente, el dato 31 es el más influyente, seguido por el 20, lo cual coincide con lo visto en el gráfico anterior.
Podríamos considerar eliminar la observación 9 para este análisis, pues no parece afectar a las componentes principales seleccionadas. Por su parte, la observación 31 parece decisiva en la elección de la tercera componente principal. El resto de datos afectan por igual a la primera y segunda componente, con algunas excepciones de datos extremos que influyen en mayor medida en una u otra componente, como por ejemplo la observación 7, con gran peso en la segunda componente pero casi nulo en la primera.
Análisis Factorial
Determinación del número óptimo de factores
Usaremos el criterio Scree plot y el análisis paralelo.
Scree Plot
Por el método del codo podemos deducir que el númeroóptimo de factores se encuentra entre 2 y 3.
Análisis paralelo
## Parallel analysis suggests that the number of factors = 2 and the number of components = NAEl análisis paralelo indica que hay al menos dos factores, siendo posible considerar un tercer factor, pues la línea de puntos del análisis paralelo cruza la línea sólida del análisis factorial en el tercer factor, o incluso el cuarto y quinto.
Lo comprobamos ahora con un test de hipótesis, que contrasta si el número de factores que se le proporcionan es o no suficiente. Le proporcionamos en primer lugar 2 factores, pues son los que parece indicar tanto el scree plot como el análisis paralelo:
factanal(datos_pca,factors=2, rotation="varimax")##
## Call:
## factanal(x = datos_pca, factors = 2, rotation = "varimax")
##
## Uniquenesses:
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI ZTLIBROP
## 0.908 0.020 0.053 0.071 0.402 0.051 0.085 0.446
## ZTEJERCI ZTPOBACT ZTENERGI
## 0.987 0.341 0.453
##
## Loadings:
## Factor1 Factor2
## ZPOBDENS 0.300
## ZTMINFAN -0.819 -0.556
## ZESPVIDA 0.830 0.508
## ZPOBURB 0.963
## ZTMEDICO 0.673 0.380
## ZPAGRICU -0.971
## ZPSERVI 0.911 -0.291
## ZTLIBROP 0.685 0.291
## ZTEJERCI 0.113
## ZTPOBACT 0.204 0.785
## ZTENERGI 0.653 0.347
##
## Factor1 Factor2
## SS loadings 5.453 1.729
## Proportion Var 0.496 0.157
## Cumulative Var 0.496 0.653
##
## Test of the hypothesis that 2 factors are sufficient.
## The chi square statistic is 63.09 on 34 degrees of freedom.
## The p-value is 0.00176El p-valor es pequeño, incluso menor que 0.01, por lo que se rechazaría la hipótesis nula de que el número de factores es suficiente. Lo ejecutamos ahora con tres factores:
factanal(datos_pca,factors=3, rotation="varimax")##
## Call:
## factanal(x = datos_pca, factors = 3, rotation = "varimax")
##
## Uniquenesses:
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI ZTLIBROP
## 0.905 0.023 0.037 0.088 0.316 0.010 0.071 0.417
## ZTEJERCI ZTPOBACT ZTENERGI
## 0.974 0.236 0.042
##
## Loadings:
## Factor1 Factor2 Factor3
## ZPOBDENS 0.295
## ZTMINFAN -0.752 -0.635
## ZESPVIDA 0.768 0.611
## ZPOBURB 0.951
## ZTMEDICO 0.619 0.405 0.370
## ZPAGRICU -0.980 -0.169
## ZPSERVI 0.943 -0.173 -0.102
## ZTLIBROP 0.637 0.339 0.249
## ZTEJERCI 0.129
## ZTPOBACT 0.777 0.390
## ZTENERGI 0.581 0.339 0.711
##
## Factor1 Factor2 Factor3
## SS loadings 5.043 1.906 0.930
## Proportion Var 0.458 0.173 0.085
## Cumulative Var 0.458 0.632 0.716
##
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 26.21 on 25 degrees of freedom.
## The p-value is 0.397En este caso el p-valor es grande, mucho mayor que 0.05, por lo que no rechazamos la hiṕotesis nula y deducimos que 3 factores son suficientes al nivel de confianza del 95%.
Todo parece indicar, por lo tanto, que el númeroóptimo de factores sería de 3, así que tomamos este número de factores para llevar a cabo el análisis factorial.
Modelo factorial
Estimamos el modelo factorial con 3 factores. Para ello, haremos una rotación usando el enfoque Varimax para poder tener una interpretación más simple. Este método maximiza la varianza por columnas de la matriz de pesos factorial.
modelo_varimax<-fa(poly_cor,nfactors = 3,rotate = "varimax",
fa="mle")Vemos la matriz de pesos factorial (rotada)
print(modelo_varimax$loadings,cut=0) ##
## Loadings:
## MR1 MR2 MR3
## ZPOBDENS 0.013 0.270 0.060
## ZTMINFAN -0.725 -0.608 -0.009
## ZESPVIDA 0.742 0.535 0.084
## ZPOBURB 0.955 0.112 0.041
## ZTMEDICO 0.607 0.531 0.179
## ZPAGRICU -0.975 -0.110 -0.023
## ZPSERVI 0.956 -0.159 -0.242
## ZTLIBROP 0.628 0.436 0.145
## ZTEJERCI 0.000 0.039 0.702
## ZTPOBACT 0.055 0.962 -0.153
## ZTENERGI 0.575 0.584 0.135
##
## MR1 MR2 MR3
## SS loadings 4.951 2.518 0.659
## Proportion Var 0.450 0.229 0.060
## Cumulative Var 0.450 0.679 0.739Veamos ahora visualmente con qué variables correlaciona cada uno de los factores
fa.diagram(modelo_varimax)
Podemos ver que el primer factor tiene una alta correlación con un gran número de variables, tanto positiva como negativa. En cambio, el segundo factor sólo se correlaciona en gran medida con dos variables, ZTPOBACT y ZTENERGI, siendo la correlación con la variable ZTPOBACT perfecta y positiva. El tercer factor tiene una correlación alta únicamente con una variable, ZTEJERCI. Hay una variable, ZPOBDENS, que presenta correlación baja con los tres factores.
En cuanto a la varianza acumulada explicada esta es de 0.739, es decir, el 74% de la varianza original queda explicada por los tres factores determinados.
Análisis de la normalidad Multivariante
outliers <- mvn(data = datos3, mvnTest = "hz", multivariateOutlierMethod = "quan")
Se detectan 11 outliers, de los cuales ya nos ocupamos anteriormente, pero ninguno de los dos test realizados a continuación encuentran evidencias al 5% de significación de falta de normalidad multivariante, es decir, al 95% de confianza podemos suponer que hay normalidad multivariante en la muestra.
Test Henze-Zirkler
hz_test <- mvn(data = datos3, mvnTest = "hz", multivariatePlot = "qq")
hz_test$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 0.9880517 0.08577614 YESTest Mardia
result <- mvn(data = datos3, mvnTest = "mardia", multivariatePlot = "qq")
result$multivariateNormality## Test Statistic p value Result
## 1 Mardia Skewness 317.293682191296 0.0983176430248101 YES
## 2 Mardia Kurtosis -0.915438362792719 0.359961546593373 YES
## 3 MVN <NA> <NA> YESAnálisis discriminante
Empezamos añadiendo una variable cualitativa con dos categorías para poder llevar a cabo la clasificación. Esta variable se ha añadido teniendo en cuenta el atributo ZESPVIDA, las instancias que queden por debajo de la mediana se han etiquetado como ‘a’ y las que queden por encima como ‘b’. Mostramos algunas instancias para ver si se ha realizado lo anterior correctamente.
labels<-c("a","a","b","b","a","b","b","a","a","a","b","a","b","b","a","a","a","b","b","b","a","a","a","a","a","b","b","b","b","b","a","b","b","a")
datos<-data.frame(datos3_nout,labels)
datos$labels<-as.factor(datos$labels)
head(datos)## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI
## 1 -0.8571914 0.9614839 -1.5498548 -0.2148886 -0.7338252 -0.60635548 0.48366865
## 2 -1.0167171 1.2133839 -0.8948879 -0.7652623 -0.8876598 0.05765957 0.05831739
## 3 -0.9918559 -0.3185794 0.4388628 1.1381136 1.2660251 -0.75993719 0.78468647
## 4 -1.0778341 -1.0203009 1.0819213 1.2802934 0.4968519 -1.05806639 1.40635370
## 5 -0.9493848 0.6119084 -0.3232805 0.4409735 -0.3107799 0.04410824 0.11066832
## 6 -1.0726547 -1.0280121 1.1295552 0.8170622 0.5160813 -1.10775460 1.62884513
## ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI labels
## 1 -0.5751040 -0.5604298 -0.7828970 0.1281385 a
## 2 -0.9696156 -0.2033670 -2.1340940 -0.5261359 a
## 3 -0.4763677 -0.4830344 -0.4189086 -0.3792220 b
## 4 1.0594990 -0.6117037 0.6244667 1.1899325 b
## 5 -0.5320795 -0.7149139 -0.4272645 -0.7252616 a
## 6 1.5770836 -0.8658606 0.9660097 -0.1690430 bExploración gráfica de los datos
Exploramos como de bien (o mal) clasifica cada una de las variables. Para ello vemos qué pares de variables separan mejor entre las dos clases.
pairs(x = datos[, c("ZPOBDENS","ZTMINFAN","ZESPVIDA","ZPOBURB", "ZTMEDICO", "ZPAGRICU", "ZPSERVI","ZTLIBROP","ZTEJERCI","ZTPOBACT","ZTENERGI")], col = c("firebrick", "green3")[datos$labels], pch = 19)
Las variables “Esperanza de vida” y “Población del sector agrícola”, por ejemplo, separan adecuadamente las dos clases.
Homogeneidad de la varianza
Contrastamos ahora la hipótesis de que la matriz de covarianzas es constante en todas las clases. Para ello usamos el test de Box M, que es una extensión del de Barttlet para escenarios multivariantes. Este test es sensible a que los datos efectivamente se distribuyan según una normal multivariante, lo cual hemos comprobado que en nuestro caso se cumple. Por este motivo se recomienda utilizar una significación (p-value) menor que 0.001.
La hipóstesis nula a contrastar es la de igualdad de matrices de covarianzas en todos los grupos.
## Loading required package: MASS## ---
## biotools version 4.2boxM(data = datos[, 1:11], grouping = datos[, 12])##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: datos[, 1:11]
## Chi-Sq (approx.) = 144.44, df = 66, p-value = 8.669e-08Puesto que el p-valor es muy pequeño, menor a 0.001, rechazamos la hipótesis nula de que las matrices de covarianzas son iguales para todos los grupos, esto es, asumimos la NO homogeneidad de varianzas. Por tanto, en este caso es conveniente usar un análisis discriminante cuadrático (QDA).
Función discriminante (cuadrática)
Se va a usar la función qda que viene en el paquete MASS
modelo_qda <- qda(formula = labels~ . ,data = datos)
modelo_qda## Call:
## qda(labels ~ ., data = datos)
##
## Prior probabilities of groups:
## a b
## 0.5 0.5
##
## Group means:
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI
## a -0.176386880 0.8279738 -0.826239 -0.7191281 -0.8282751 0.747183 -0.5641197
## b 0.002955307 -0.8279738 0.826239 0.7191281 0.8282751 -0.747183 0.5641197
## ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## a -0.6735376 -0.3107902 -0.5583716 -0.7157038
## b 0.6735376 -0.2685980 0.5583716 0.3776178La salida de este objeto, nos muestra las probabilidades a priori de cada grupo, en este caso 0.5 y las medias de cada regresor por grupo.
Una vez construido el clasificador podemos clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict. Por ejemplo, veamos como se clasifica un nuevo registros con los siguientes valores ZPOBDENS=0.3,ZTMINFAN=-0.2,ZESPVIDA=0.1,ZPOBURB=-0.1,ZTMEDICO=0.5,ZPAGRICU=-0.5, ZPSERVI=0.5, ZTLIBROP=0.3, ZTEJERCI=-0.8,ZTPOBACT=0.2,ZTENERGI=0.1
nuevas_observaciones <- data.frame(ZPOBDENS=0.3,ZTMINFAN=-0.2,ZESPVIDA=0.1,ZPOBURB=-0.1,ZTMEDICO=0.5,ZPAGRICU=-0.5, ZPSERVI=0.5, ZTLIBROP=0.3, ZTEJERCI=-0.8,ZTPOBACT=0.2,ZTENERGI=0.1)
predict(object = modelo_qda, newdata = nuevas_observaciones)## $class
## [1] a
## Levels: a b
##
## $posterior
## a b
## 1 0.917501 0.08249898Según la función discriminante, la probabilidad a posteriori de que el nuevo registro tenga la etiqueta a es del 91,7% mientras la de que esté en la etiqueta b es inferior al 1%. Por tanto este dato sería clasificado con la etiqueta a. En efecto, el valor de ZESPVIDA para este nuevo registro (0.1) está por debajo de la mediana (0.2781) y este hecho tiene mucha influencia.
Matriz de confusión
Construimos ahora una matriz de confusión usando previamente validación cruzada en el modelo de clasificación
library('biotools')
pred <- predict(modelo_qda, dimen = 1)
confusionmatrix(datos$labels, pred$class)## new a new b
## a 17 0
## b 0 17Podemos comprobar que todos los ejemplos de entrenamiento se clasifican correctamente. Veamos ahora el porcentaje de errores de clasificación:
trainig_error <- mean(datos$labels != pred$class) * 100
paste("trainig_error=", trainig_error, "%")## [1] "trainig_error= 0 %"El porcentaje de aciertos es del 100%. Este hecho nos dice que el entrenamiento ha sido perfecto, lo que nos lleva a pensar en la posibilidad de un sobreajuste, aunque con tampoco datos y con la facilidad de la base de datos podemos descartarlo.
Visualización de las clasificaciones
Representamos los límites de clasificación del modelo para cada par de predictores. Cada color representa una región de clasificación acorde al modelo, se muestra el centroide de cada región y el valor real de las observaciones.
par(mar=c(1,1,1,1))
partimat(formula =labels ~ datos$ZPOBDENS + datos$ZTMINFAN + datos$ZESPVIDA + datos$ZPOBURB , data = datos,
method = "qda", prec = 400,
image.colors = c("darkgoldenrod1", "skyblue2"),
col.mean = "firebrick")
Observamos que hay varios pares de variables que clasifican perfectamente todas las observaciones, por lo que quizás no sería necesario considerar todas las variables.
Debido al número de datos tan reducido del que disponemos, quizás llevar a cabo un análisis cuadrático no tenga mucho sentido, a pesar de la no homegenidad de la varianza ante la que nos encontramos. Por tanto, realizaremos un análisis discriminante lineal para comparar los resultados.
Función discriminante (lineal)
Ahora vamos a aplicar el análisis discriminante lineal, para ello hay que tener en cuenta que todas las variables deben distribuirse normalmente ya que de lo contrario los resultados estarán sesgados. Además también se debe de tener en cuenta los outliers y que la escala de todas las variables sean comparables, para ello se estandarizan. El método LDA asume que los predictores están normalmente distribuidos, es decir, provienen de la distribución gaussiana.
Creamos el modelo y lo entrenamos
modelo_lda<-lda(formula = labels~ . ,data = datos)
modelo_lda## Call:
## lda(labels ~ ., data = datos)
##
## Prior probabilities of groups:
## a b
## 0.5 0.5
##
## Group means:
## ZPOBDENS ZTMINFAN ZESPVIDA ZPOBURB ZTMEDICO ZPAGRICU ZPSERVI
## a -0.176386880 0.8279738 -0.826239 -0.7191281 -0.8282751 0.747183 -0.5641197
## b 0.002955307 -0.8279738 0.826239 0.7191281 0.8282751 -0.747183 0.5641197
## ZTLIBROP ZTEJERCI ZTPOBACT ZTENERGI
## a -0.6735376 -0.3107902 -0.5583716 -0.7157038
## b 0.6735376 -0.2685980 0.5583716 0.3776178
##
## Coefficients of linear discriminants:
## LD1
## ZPOBDENS -0.02696809
## ZTMINFAN 0.63418453
## ZESPVIDA 1.16461874
## ZPOBURB -0.19204910
## ZTMEDICO 0.91738184
## ZPAGRICU -1.07071939
## ZPSERVI -0.04748110
## ZTLIBROP 0.17210741
## ZTEJERCI -0.18187509
## ZTPOBACT 0.60599703
## ZTENERGI 0.03187244Una vez construido el clasificador podemos clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict. Por ejemplo, veamos como se clasifica un nuevo registros con los siguientes valores ZPOBDENS=0.3,ZTMINFAN=-0.2,ZESPVIDA=0.1,ZPOBURB=-0.1,ZTMEDICO=0.5,ZPAGRICU=-0.5, ZPSERVI=0.5, ZTLIBROP=0.3, ZTEJERCI=-0.8,ZTPOBACT=0.2,ZTENERGI=0.1
nuevas_observaciones <- data.frame(ZPOBDENS=0.3,ZTMINFAN=-0.2,ZESPVIDA=0.1,ZPOBURB=-0.1,ZTMEDICO=0.5,ZPAGRICU=-0.5, ZPSERVI=0.5, ZTLIBROP=0.3, ZTEJERCI=-0.8,ZTPOBACT=0.2,ZTENERGI=0.1)
predict(object = modelo_lda, newdata = nuevas_observaciones)## $class
## [1] b
## Levels: a b
##
## $posterior
## a b
## 1 0.003296763 0.9967032
##
## $x
## LD1
## 1 1.24293Observamos ahora que la etiqueta predicha para este nuevo registro es la b, entrando en conflicto con la predicción por parte del predictor cuadrático. Aún así tampoco podemos saber con certeza que predicción es la correcta y cual la incorrecta, según el criterio usado en la elección de las etiquetas, es lógico pensar que la etiqueta correcta sería la a y la incorrecta la b. Pero puede ocurrir que existan otro tipo de relaciones entre atributos que no hayamos tenido en cuenta y que influyan, provocando este conflicto entre las etiquetas predichas por parte de ambos predictores.
Matriz de confusión
Mostramos a continuación la matriz de confusión usando validación cruzada en los ejemplos de entrenamiento
pred <- predict(modelo_lda, dimen = 1)
confusionmatrix(datos$labels, pred$class)## new a new b
## a 17 0
## b 1 16Podemos comprobar que únicamente un dato de entrenamiento se clasifican incorrectamente. Veamos ahora el porcentaje de errores de clasificación:
trainig_error <- mean(datos$labels != pred$class) * 100
paste("trainig_error=", trainig_error, "%")## [1] "trainig_error= 2.94117647058824 %"En consecuencia, el porcentaje de aciertos es del 97.1%. Podemos notar que la clasificación en el entrenamiento no ha sido perfecta pues es lógico pensar que la clasificación anterior pudiera haber sido errónea. Sin embargo, al ser la tasa de acierto tan alta tampoco resulta fácil afirmarlo con behemencia. Por otro lado, al no tener un acierto del 100%, podemos llegar a pensar que la capacidad de generalizar puede ser mayor que el caso anterior, pues no habría tanto sobreajuste. Aunque como antes se explicó, podemos llegar a descartar la presencia de sobreajuste.
Visualización de las clasificaciones
Como hicimos en el análisis discriminante con la función cuadrática, aquí también se representará los límites de clasificación del modelo para cada par de predictores. Cada color representa una región de clasificación acorde al modelo, se muestra el centroide de cada región y el valor real de las observaciones.
par(mar=c(1,1,1,1))
partimat(formula =labels ~ datos$ZPOBDENS + datos$ZTMINFAN + datos$ZESPVIDA + datos$ZPOBURB , data = datos,
method = "lda", prec = 400,
image.colors = c("darkgoldenrod1", "skyblue2"),
col.mean = "firebrick")
Como conclusión general, podemos decir que en el entrenamiento QDA lo hace perfecto mientras que LDA es casi perfecto, por lo que pueden existir ciertos indicios de sobreajuste siendo esto más notorio en QDA que en LDA, aunque debido a las pocas instancias de la base de datos y a este ejemplo que nos hemos inventado podríamos descartar esta última idea. En la predicción ambos predictores entran en conflicto con un mismo registro, siendo la predicción del QDA más razonable que la del LDA aunque no se puede comprobar a ciencia cierta cual ha acertado y cual ha fallado, pues nos hemos inventado las etiquetas de acuerdo a un criterio sin saber cuán de lógico es.
Teniendo en cuenta que los límites de decisión del LDA y del QDA son lineales y cuadráticos respectivamente, podemos deducir que para grandes bases de datos, el QDA es más factible pues tiene mayor flexibilidad que el LDA, presentando este menos sesgo pero mayor ajuste a la varianza de los datos. Mientras que para bases de datos más reducidas, como por ejemplo esta, LDA debería conseguir mejores clasificaciones de nuevos datos, esto es, generaliza mejor en datos nuevos que QDA, pues al entrenar el modelo con pocos ejemplos de entrenamiento se nos presenta un escenario en el que evitar la varianza es crucial.